Data analysis for:
Why are telomeres the length that they are? Insight from a phylogenetic comparative analysisWorkflow of the analyses produced by Dylan Padilla, School of Life Sciences, Arizona State University, Tempe, AZ 85287, USA.
Library
library(AICcmodavg)
library(ape)
library(caper)
library(car)
library(coda)
library(extrafont)
library(geiger)
library(graph)
library(kableExtra)
library(MuMIn)
library(nlme)
library(pbapply)
library(phylopath)
library(phytools)
library(plotrix)
library(rphylopic)
library(scales)
library(sensiPhy)
library(shape)
library(xtable)R.version> _
> platform aarch64-apple-darwin20
> arch aarch64
> os darwin20
> system aarch64, darwin20
> status
> major 4
> minor 4.1
> year 2024
> month 06
> day 14
> svn rev 86737
> language R
> version.string R version 4.4.1 (2024-06-14)
> nickname Race for Your LifesessionInfo()> R version 4.4.1 (2024-06-14)
> Platform: aarch64-apple-darwin20
> Running under: macOS Sonoma 14.6.1
>
> Matrix products: default
> BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
>
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
>
> time zone: America/Phoenix
> tzcode source: internal
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] rmarkdown_2.27 xtable_1.8-4 shape_1.4.6.1 sensiPhy_0.8.5 ggplot2_3.5.1 phylolm_2.6.2
> [7] scales_1.3.0 rphylopic_1.4.0.9000 plotrix_3.8-4 phylopath_1.3.0 pbapply_1.7-2 nlme_3.1-165
> [13] MuMIn_1.48.4 kableExtra_1.4.0 graph_1.82.0 BiocGenerics_0.50.0 geiger_2.0.11 phytools_2.3-0
> [19] maps_3.4.2 extrafont_0.19 coda_0.19-4.1 car_3.1-2 carData_3.0-5 caper_1.0.3
> [25] mvtnorm_1.2-5 MASS_7.3-60.2 ape_5.8 AICcmodavg_2.3-3
>
> loaded via a namespace (and not attached):
> [1] mnormt_2.1.1 gridExtra_2.3 phangorn_2.11.1 rlang_1.1.4 magrittr_2.0.3
> [6] compiler_4.4.1 png_0.1-8 systemfonts_1.1.0 vctrs_0.6.5 combinat_0.0-8
> [11] quadprog_1.5-8 stringr_1.5.1 pkgconfig_2.0.3 fastmap_1.2.0 ggraph_2.2.1
> [16] labeling_0.4.3 utf8_1.2.4 subplex_1.9 deSolve_1.40 purrr_1.0.2
> [21] xfun_0.45 cachem_1.1.0 clusterGeneration_1.3.8 jsonlite_1.8.8 highr_0.11
> [26] tweenr_2.0.3 jpeg_0.1-10 VGAM_1.1-11 parallel_4.4.1 R6_2.5.1
> [31] bslib_0.7.0 stringi_1.8.4 parallelly_1.38.0 extrafontdb_1.0 jquerylib_0.1.4
> [36] numDeriv_2016.8-1.1 Rcpp_1.0.12 iterators_1.0.14 knitr_1.47 future.apply_1.11.2
> [41] optimParallel_1.0-2 ggm_2.5.1 base64enc_0.1-3 Matrix_1.7-0 splines_4.4.1
> [46] igraph_2.0.3 tidyselect_1.2.1 viridis_0.6.5 yaml_2.3.8 rstudioapi_0.16.0
> [51] abind_1.4-5 doParallel_1.0.17 codetools_0.2-20 curl_5.2.1 listenv_0.9.1
> [56] lattice_0.22-6 tibble_3.2.1 withr_3.0.0 evaluate_0.24.0 future_1.34.0
> [61] survival_3.6-4 polyclip_1.10-7 xml2_1.3.6 BiocManager_1.30.23 pillar_1.9.0
> [66] foreach_1.5.2 stats4_4.4.1 generics_0.1.3 grImport2_0.3-1 munsell_0.5.1
> [71] globals_0.16.3 glue_1.7.0 unmarked_1.4.1 scatterplot3d_0.3-44 tools_4.4.1
> [76] graphlayouts_1.1.1 XML_3.99-0.17 tidygraph_1.3.1 fastmatch_1.1-4 grid_4.4.1
> [81] tidyr_1.3.1 Rttf2pt1_1.3.12 colorspace_2.1-0 ggforce_0.4.2 cli_3.6.3
> [86] DEoptim_2.2-8 rsvg_2.6.0 fansi_1.0.6 expm_0.999-9 viridisLite_0.4.2
> [91] svglite_2.1.3 dplyr_1.1.4 gtable_0.3.5 sass_0.4.9 digest_0.6.36
> [96] ggrepel_0.9.5 farver_2.1.2 memoise_2.0.1 htmltools_0.5.8.1 lifecycle_1.0.4
> [101] httr_1.4.7set.seed(80)
## Dataset
data <- read.csv('../TRF-anaysis-no-domesticated.csv')
str(data)> 'data.frame': 131 obs. of 27 variables:
> $ Common.name : chr "European seabass" "Skipjack trevally" "Green swordtail" "Southern platyfish" ...
> $ Scientific_name : chr "Dicentrarchus_labrax" "Pseudocaranx_wrighti" "Xiphophorus_hellerii" "Xiphophorus_maculatus" ...
> $ Domain : chr "Eukaryota" "Eukaryota" "Eukaryota" "Eukaryota" ...
> $ Kingdom : chr "Metazoa" "Metazoa" "Metazoa" "Metazoa" ...
> $ Phylum : chr "Chordata" "Chordata" "Chordata" "Chordata" ...
> $ Class : chr "Fish" "Fish" "Fish" "Fish" ...
> $ Order : chr "Moroniformes" "Carangiformes" "Cyprinodontiformes" "Cyprinodontiformes" ...
> $ Family : chr "Moronidae" "Carangidae" "Poeciliidae" "Poeciliidae" ...
> $ Genus : chr "Dicentrarchus" "Pseudocaranx" "Xiphophorus" "Xiphophorus" ...
> $ Species : chr "labrax" "wrighti" "hellerii" "maculatus" ...
> $ Endo_ectotherm : chr "ecto" "ecto" "ecto" "ecto" ...
> $ Adult_mass_grams : num 6600 3000 3.7 1.8 1600 ...
> $ log_mass : num 3.82 3.477 0.568 0.255 3.204 ...
> $ Lifespan_years : num 15 13 5 5 41 47 18 10 7 20 ...
> $ Loglifespan : num 1.176 1.114 0.699 0.699 1.613 ...
> $ Average_Telomere_Length_kb: num 3.44 4.06 4.25 4.25 4.6 4.62 5.26 5.48 5.56 6.14 ...
> $ Telomerase_activity : int NA NA NA NA NA NA NA NA NA NA ...
> $ Tissue.type.for.TA : chr NA NA NA NA ...
> $ Source.for.TA : chr NA NA NA NA ...
> $ Tissue.type.for.TL : chr "Blood" "Muscle" "Multiple somatic tissues" "Multiple somatic tissues" ...
> $ Methodology.for.TL : chr "TRF (denatured)" "TRF (denatured)" "TRF (denatured)" "TRF (denatured)" ...
> $ TRF_type : chr "A" "A" "A" "A" ...
> $ Source.for.TL : chr "Horn et al. 2008" "Izzo 2010" "Downs et al. 2012" "Downs et al. 2012" ...
> $ Source.for.Lifespan : chr "AnAge" "Ratnasingham and Herbert 2007" "AnAge" "National park aquarium" ...
> $ Source.for.Mass : chr "AnAge" "Marinewise" "Meyer et al. 2006" "McKenzie et al. 1983" ...
> $ Captivity : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Notes : chr "" "" "" "" ...head(data)> Common.name Scientific_name Domain Kingdom Phylum Class Order Family Genus Species
> 1 European seabass Dicentrarchus_labrax Eukaryota Metazoa Chordata Fish Moroniformes Moronidae Dicentrarchus labrax
> 2 Skipjack trevally Pseudocaranx_wrighti Eukaryota Metazoa Chordata Fish Carangiformes Carangidae Pseudocaranx wrighti
> 3 Green swordtail Xiphophorus_hellerii Eukaryota Metazoa Chordata Fish Cyprinodontiformes Poeciliidae Xiphophorus hellerii
> 4 Southern platyfish Xiphophorus_maculatus Eukaryota Metazoa Chordata Fish Cyprinodontiformes Poeciliidae Xiphophorus maculatus
> 5 Goldfish Carassius_auratus Eukaryota Metazoa Chordata Fish Cypriniformes Cyprinidae Carassius auratus
> 6 Common carp Cyprinus_carpio Eukaryota Metazoa Chordata Fish Cypriniformes Cyprinidae Cyprinus carpio
> Endo_ectotherm Adult_mass_grams log_mass Lifespan_years Loglifespan Average_Telomere_Length_kb Telomerase_activity Tissue.type.for.TA
> 1 ecto 6600.0 3.8195439 15 1.176091 3.44 NA <NA>
> 2 ecto 3000.0 3.4771213 13 1.113943 4.06 NA <NA>
> 3 ecto 3.7 0.5682017 5 0.698970 4.25 NA <NA>
> 4 ecto 1.8 0.2552725 5 0.698970 4.25 NA <NA>
> 5 ecto 1600.0 3.2041200 41 1.612784 4.60 NA <NA>
> 6 ecto 10600.0 4.0253059 47 1.672098 4.62 NA <NA>
> Source.for.TA Tissue.type.for.TL Methodology.for.TL TRF_type Source.for.TL Source.for.Lifespan Source.for.Mass
> 1 <NA> Blood TRF (denatured) A Horn et al. 2008 AnAge AnAge
> 2 <NA> Muscle TRF (denatured) A Izzo 2010 Ratnasingham and Herbert 2007 Marinewise
> 3 <NA> Multiple somatic tissues TRF (denatured) A Downs et al. 2012 AnAge Meyer et al. 2006
> 4 <NA> Multiple somatic tissues TRF (denatured) A Downs et al. 2012 National park aquarium McKenzie et al. 1983
> 5 <NA> Muscle TRF (denatured) A Izzo 2010 AnAge Fishbase
> 6 <NA> Muscle TRF (denatured) A Izzo 2010 AnAge AnAge
> Captivity Notes
> 1 0
> 2 0
> 3 0
> 4 0
> 5 0
> 6 0unique(data$Class)> [1] "Fish" "Reptilia" "Mammalia" "Aves"dat <- data
names(dat)> [1] "Common.name" "Scientific_name" "Domain" "Kingdom"
> [5] "Phylum" "Class" "Order" "Family"
> [9] "Genus" "Species" "Endo_ectotherm" "Adult_mass_grams"
> [13] "log_mass" "Lifespan_years" "Loglifespan" "Average_Telomere_Length_kb"
> [17] "Telomerase_activity" "Tissue.type.for.TA" "Source.for.TA" "Tissue.type.for.TL"
> [21] "Methodology.for.TL" "TRF_type" "Source.for.TL" "Source.for.Lifespan"
> [25] "Source.for.Mass" "Captivity" "Notes"unique(is.na(dat$log_mass))> [1] FALSE## Trees
##write.table(data$Scientific_name, 'spp.tree.08-03-24.txt', row.names = FALSE, col.names = FALSE)
full_data_tree <- read.tree("spp.tree.08-03-24.nwk")
is.ultrametric(full_data_tree)> [1] FALSEfull_data_tree <- force.ultrametric(full_data_tree)> ***************************************************************
> * Note: *
> * force.ultrametric does not include a formal method to *
> * ultrametricize a tree & should only be used to coerce *
> * a phylogeny that fails is.ultrametric due to rounding -- *
> * not as a substitute for formal rate-smoothing methods. *
> ***************************************************************is.ultrametric(full_data_tree)> [1] TRUEfull_data_tree>
> Phylogenetic tree with 138 tips and 137 internal nodes.
>
> Tip labels:
> Urolophus_paucimaculatus, Callorhinchus_milii, Oncorhynchus_mykiss, Lutjanus_argentimaculatus, Acanthopagrus_schlegelii, Dicentrarchus_labrax, ...
> Node labels:
> , '14', '12', '48', '57', '49', ...
>
> Rooted; includes branch lengths.rownames(dat) <- dat$Scientific_name
## Full_tree
check <- name.check(full_data_tree, dat)
rm_phy <- check$tree_not_data
rm_dat <- check$data_not_tree
pruned_data_tree <- drop.tip(full_data_tree, rm_phy)
pruned_dat <- subset(dat, subset = dat$Scientific_name %in% pruned_data_tree$tip, select = names(dat))
str(pruned_dat)> 'data.frame': 123 obs. of 27 variables:
> $ Common.name : chr "European seabass" "Skipjack trevally" "Green swordtail" "Southern platyfish" ...
> $ Scientific_name : chr "Dicentrarchus_labrax" "Pseudocaranx_wrighti" "Xiphophorus_hellerii" "Xiphophorus_maculatus" ...
> $ Domain : chr "Eukaryota" "Eukaryota" "Eukaryota" "Eukaryota" ...
> $ Kingdom : chr "Metazoa" "Metazoa" "Metazoa" "Metazoa" ...
> $ Phylum : chr "Chordata" "Chordata" "Chordata" "Chordata" ...
> $ Class : chr "Fish" "Fish" "Fish" "Fish" ...
> $ Order : chr "Moroniformes" "Carangiformes" "Cyprinodontiformes" "Cyprinodontiformes" ...
> $ Family : chr "Moronidae" "Carangidae" "Poeciliidae" "Poeciliidae" ...
> $ Genus : chr "Dicentrarchus" "Pseudocaranx" "Xiphophorus" "Xiphophorus" ...
> $ Species : chr "labrax" "wrighti" "hellerii" "maculatus" ...
> $ Endo_ectotherm : chr "ecto" "ecto" "ecto" "ecto" ...
> $ Adult_mass_grams : num 6600 3000 3.7 1.8 1600 ...
> $ log_mass : num 3.82 3.477 0.568 0.255 3.204 ...
> $ Lifespan_years : num 15 13 5 5 41 47 18 10 7 20 ...
> $ Loglifespan : num 1.176 1.114 0.699 0.699 1.613 ...
> $ Average_Telomere_Length_kb: num 3.44 4.06 4.25 4.25 4.6 4.62 5.26 5.48 5.56 6.14 ...
> $ Telomerase_activity : int NA NA NA NA NA NA NA NA NA NA ...
> $ Tissue.type.for.TA : chr NA NA NA NA ...
> $ Source.for.TA : chr NA NA NA NA ...
> $ Tissue.type.for.TL : chr "Blood" "Muscle" "Multiple somatic tissues" "Multiple somatic tissues" ...
> $ Methodology.for.TL : chr "TRF (denatured)" "TRF (denatured)" "TRF (denatured)" "TRF (denatured)" ...
> $ TRF_type : chr "A" "A" "A" "A" ...
> $ Source.for.TL : chr "Horn et al. 2008" "Izzo 2010" "Downs et al. 2012" "Downs et al. 2012" ...
> $ Source.for.Lifespan : chr "AnAge" "Ratnasingham and Herbert 2007" "AnAge" "National park aquarium" ...
> $ Source.for.Mass : chr "AnAge" "Marinewise" "Meyer et al. 2006" "McKenzie et al. 1983" ...
> $ Captivity : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Notes : chr "" "" "" "" ...head(pruned_dat)> Common.name Scientific_name Domain Kingdom Phylum Class Order Family
> Dicentrarchus_labrax European seabass Dicentrarchus_labrax Eukaryota Metazoa Chordata Fish Moroniformes Moronidae
> Pseudocaranx_wrighti Skipjack trevally Pseudocaranx_wrighti Eukaryota Metazoa Chordata Fish Carangiformes Carangidae
> Xiphophorus_hellerii Green swordtail Xiphophorus_hellerii Eukaryota Metazoa Chordata Fish Cyprinodontiformes Poeciliidae
> Xiphophorus_maculatus Southern platyfish Xiphophorus_maculatus Eukaryota Metazoa Chordata Fish Cyprinodontiformes Poeciliidae
> Carassius_auratus Goldfish Carassius_auratus Eukaryota Metazoa Chordata Fish Cypriniformes Cyprinidae
> Cyprinus_carpio Common carp Cyprinus_carpio Eukaryota Metazoa Chordata Fish Cypriniformes Cyprinidae
> Genus Species Endo_ectotherm Adult_mass_grams log_mass Lifespan_years Loglifespan
> Dicentrarchus_labrax Dicentrarchus labrax ecto 6600.0 3.8195439 15 1.176091
> Pseudocaranx_wrighti Pseudocaranx wrighti ecto 3000.0 3.4771213 13 1.113943
> Xiphophorus_hellerii Xiphophorus hellerii ecto 3.7 0.5682017 5 0.698970
> Xiphophorus_maculatus Xiphophorus maculatus ecto 1.8 0.2552725 5 0.698970
> Carassius_auratus Carassius auratus ecto 1600.0 3.2041200 41 1.612784
> Cyprinus_carpio Cyprinus carpio ecto 10600.0 4.0253059 47 1.672098
> Average_Telomere_Length_kb Telomerase_activity Tissue.type.for.TA Source.for.TA Tissue.type.for.TL
> Dicentrarchus_labrax 3.44 NA <NA> <NA> Blood
> Pseudocaranx_wrighti 4.06 NA <NA> <NA> Muscle
> Xiphophorus_hellerii 4.25 NA <NA> <NA> Multiple somatic tissues
> Xiphophorus_maculatus 4.25 NA <NA> <NA> Multiple somatic tissues
> Carassius_auratus 4.60 NA <NA> <NA> Muscle
> Cyprinus_carpio 4.62 NA <NA> <NA> Muscle
> Methodology.for.TL TRF_type Source.for.TL Source.for.Lifespan Source.for.Mass Captivity Notes
> Dicentrarchus_labrax TRF (denatured) A Horn et al. 2008 AnAge AnAge 0
> Pseudocaranx_wrighti TRF (denatured) A Izzo 2010 Ratnasingham and Herbert 2007 Marinewise 0
> Xiphophorus_hellerii TRF (denatured) A Downs et al. 2012 AnAge Meyer et al. 2006 0
> Xiphophorus_maculatus TRF (denatured) A Downs et al. 2012 National park aquarium McKenzie et al. 1983 0
> Carassius_auratus TRF (denatured) A Izzo 2010 AnAge Fishbase 0
> Cyprinus_carpio TRF (denatured) A Izzo 2010 AnAge AnAge 0pruned_data_tree>
> Phylogenetic tree with 123 tips and 122 internal nodes.
>
> Tip labels:
> Urolophus_paucimaculatus, Callorhinchus_milii, Oncorhynchus_mykiss, Lutjanus_argentimaculatus, Acanthopagrus_schlegelii, Dicentrarchus_labrax, ...
> Node labels:
> , '14', '12', '48', '57', '49', ...
>
> Rooted; includes branch lengths.name.check(pruned_data_tree, pruned_dat)> [1] "OK"## Checking the distribution of lifespan
plot(NA, xlim = c(0, 200), ylim = c(0, 60), type = 'n', las = 1, xlab = '', ylab = '', axes = FALSE)
grid()
par(new = TRUE)
hist(pruned_dat$Lifespan_years, main = "Raw variable", las = 1, xlab = 'Lifespan (years)')
box()
plot(NA, xlim = c(0, 200), ylim = c(0, 60), type = 'n', las = 1, xlab = '', ylab = '', axes = FALSE)
grid()
par(new = TRUE)
hist(log1p(pruned_dat$Lifespan_years), main = "Log-transformed", xlab = 'Log lifespan (years)', las = 1)
box()
## Log-transforming for a better fit
pruned_dat$log.lifespan <- log1p(pruned_dat$Lifespan_years)
pruned_dat$log.mass <- log1p(pruned_dat$Adult_mass_grams)Implementing phylogenetic path analysis as suggested by the editor
## Defining models
path.mod <- pruned_dat
names(path.mod)> [1] "Common.name" "Scientific_name" "Domain" "Kingdom"
> [5] "Phylum" "Class" "Order" "Family"
> [9] "Genus" "Species" "Endo_ectotherm" "Adult_mass_grams"
> [13] "log_mass" "Lifespan_years" "Loglifespan" "Average_Telomere_Length_kb"
> [17] "Telomerase_activity" "Tissue.type.for.TA" "Source.for.TA" "Tissue.type.for.TL"
> [21] "Methodology.for.TL" "TRF_type" "Source.for.TL" "Source.for.Lifespan"
> [25] "Source.for.Mass" "Captivity" "Notes" "log.lifespan"
> [29] "log.mass"path.mod$Average_Telomere_Length_kb <- log(path.mod$Average_Telomere_Length_kb)
names(path.mod)[c(16, 11, 29, 28)] <- c('TL', 'TB', 'BM', 'LS')
models <- define_model_set(
"null" = c(TL ~ TL, TB ~ LS, TB ~ BM, LS ~ BM),
"(a)" = c(TL ~ LS, TB ~ LS, TB ~ BM, BM ~ LS),
"(b)" = c(TL ~ BM, TB ~ LS, TB ~ BM, BM ~ LS),
"(c)" = c(TL ~ TB, TB ~ LS, TB ~ BM, BM ~ LS),
"(d)" = c(TL ~ BM + LS + TB),
"(e)" = c(TL ~ LS),
"(f)" = c(TL ~ BM),
"(g)" = c(TL ~ TB))
#png("figureS1-revised.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figureS1-revised.pdf")
plot_model_set(models, nrow = 4) + ggplot2::scale_y_continuous(expand = c(0.1, 0)) + theme(plot.margin=unit(c(0.5, 3, 0.5, 3), "cm"))
#dev.off()
mo <- phylo_path(models, path.mod, pruned_data_tree)
s <- summary(mo)
theme_set(theme_bw())
#png("figureS2-revised.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figureS2-revised.pdf")
bar.pl <- tapply(s$w, list(s$model), print)> [1] 0.06045404
> [1] 0.03543754
> [1] 0.9037687
> [1] 3.83462e-08
> [1] 3.006297e-11
> [1] 9.501932e-12
> [1] 9.928878e-12
> [1] 0.0003397325p.val <- s[names(bar.pl), ]
barplot(NA, beside = TRUE, space = FALSE, horiz = TRUE, xlim = c(0, 1), axes = FALSE)
grid()
par(new = TRUE)
barplot(bar.pl, beside = TRUE, space = FALSE, horiz = TRUE, xlim = c(0, 1), yaxt = 'n', ylab = 'Candidate models', xlab = 'Weight of evidence (w)', col = c(rep('gray', 2), alpha('red', 0.3), rep('gray', 5)))
axis(2, at = (0:7)+0.5, labels = names(bar.pl), las = 1)
box()
text(x = round(bar.pl, 3)+0.05, y = (0:7)+0.5, labels = round(p.val$p, 3))
legend('topright', legend = 'within 2 CICc', fil = alpha('red', 0.3), bty = 'n')
#dev.off()
pagel.model0 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(pagel.model0)> Generalized least squares fit by maximum likelihood
> Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass
> Data: pruned_dat
> AIC BIC logLik
> 235.1572 252.0303 -111.5786
>
> Correlation Structure: corPagel
> Formula: ~Scientific_name
> Parameter estimate(s):
> lambda
> 0.6254488
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.629715 0.3916907 9.266787 0.0000
> Endo_ectothermendo -0.819368 0.2207192 -3.712263 0.0003
> log.lifespan -0.134538 0.1130636 -1.189932 0.2364
> log.mass -0.033764 0.0258447 -1.306413 0.1939
>
> Correlation:
> (Intr) End_ct lg.lfs
> Endo_ectothermendo -0.041
> log.lifespan -0.502 -0.101
> log.mass 0.075 0.013 -0.678
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -2.34917752 -0.02875553 0.50192730 0.92152006 2.39155282
>
> Residual standard error: 0.8117616
> Degrees of freedom: 123 total; 119 residualbrown.model <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass, correlation = corBrownian(phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(brown.model)> Generalized least squares fit by maximum likelihood
> Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass
> Data: pruned_dat
> AIC BIC logLik
> 272.3855 286.4464 -131.1927
>
> Correlation Structure: corBrownian
> Formula: ~Scientific_name
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.901248 0.8933978 4.366753 0.0000
> Endo_ectothermendo -1.172484 0.2869553 -4.085949 0.0001
> log.lifespan -0.142480 0.1268266 -1.123421 0.2635
> log.mass -0.048054 0.0328901 -1.461046 0.1466
>
> Correlation:
> (Intr) End_ct lg.lfs
> Endo_ectothermendo -0.025
> log.lifespan -0.268 -0.136
> log.mass -0.086 0.108 -0.499
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.11576421 0.04130495 0.29665511 0.48025508 1.15040877
>
> Residual standard error: 1.926538
> Degrees of freedom: 123 total; 119 residualVerTree <- pruned_data_tree
VerTree$edge.length <- pruned_data_tree$edge.length * 100
OU.model1 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass, correlation = corMartins(1, phy = VerTree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(OU.model1)> Generalized least squares fit by maximum likelihood
> Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass
> Data: pruned_dat
> AIC BIC logLik
> 283.5872 300.4603 -135.7936
>
> Correlation Structure: corMartins
> Formula: ~Scientific_name
> Parameter estimate(s):
> alpha
> 1
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.0127181 0.27461330 10.970765 0.0000
> Endo_ectothermendo -0.0099686 0.15910482 -0.062654 0.9501
> log.lifespan -0.1146562 0.11583847 -0.989794 0.3243
> log.mass -0.0007418 0.02556906 -0.029011 0.9769
>
> Correlation:
> (Intr) End_ct lg.lfs
> Endo_ectothermendo -0.108
> log.lifespan -0.748 -0.291
> log.mass 0.131 0.098 -0.650
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -3.2551171 -0.5625105 -0.0212407 0.5130302 1.9228926
>
> Residual standard error: 0.7298431
> Degrees of freedom: 123 total; 119 residual## Model selection procedure based on AIC criterion
output.models <- model.sel(brown.model, pagel.model0, OU.model1)
sel.table <- round(as.data.frame(output.models)[6:10], 3)
names(sel.table)[1] <- "K"
sel.table$Model <- rownames(sel.table)
sel.table <- sel.table[ , c(6, 1, 2, 3, 4, 5)]
sel.table> Model K logLik AICc delta weight
> pagel.model0 pagel.model0 6 -111.579 235.881 0.000 1
> brown.model brown.model 5 -131.193 272.898 37.017 0
> OU.model1 OU.model1 6 -135.794 284.311 48.430 0anova(pagel.model0)> Denom. DF: 119
> numDF F-value p-value
> (Intercept) 1 84.64587 <.0001
> Endo_ectotherm 1 16.64769 0.0001
> log.lifespan 1 7.96289 0.0056
> log.mass 1 1.70672 0.1939pagel.mod.red <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(pagel.mod.red)> Generalized least squares fit by maximum likelihood
> Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan
> Data: pruned_dat
> AIC BIC logLik
> 234.9012 248.9622 -112.4506
>
> Correlation Structure: corPagel
> Formula: ~Scientific_name
> Parameter estimate(s):
> lambda
> 0.6138179
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.662657 0.3875853 9.449939 0.0000
> Endo_ectothermendo -0.809681 0.2203418 -3.674661 0.0004
> log.lifespan -0.233895 0.0830949 -2.814794 0.0057
>
> Correlation:
> (Intr) End_ct
> Endo_ectothermendo -0.043
> log.lifespan -0.620 -0.126
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -2.293252813 -0.003891462 0.546445007 0.947511030 2.262724689
>
> Residual standard error: 0.8098694
> Degrees of freedom: 123 total; 120 residualanova(pagel.model0, pagel.mod.red) # no difference between the models> Model df AIC BIC logLik Test L.Ratio p-value
> pagel.model0 1 6 235.1572 252.0303 -111.5786
> pagel.mod.red 2 5 234.9013 248.9622 -112.4506 1 vs 2 1.744052 0.1866mod.plot <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm*log.lifespan, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML") ## for plotting purposes. Notice the slopes are not significant, but the main effects are!
summary(mod.plot)> Generalized least squares fit by maximum likelihood
> Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.lifespan
> Data: pruned_dat
> AIC BIC logLik
> 235.0328 251.9059 -111.5164
>
> Correlation Structure: corPagel
> Formula: ~Scientific_name
> Parameter estimate(s):
> lambda
> 0.5521672
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.222911 0.4660986 6.914656 0.0000
> Endo_ectothermendo -0.103630 0.5207431 -0.199003 0.8426
> log.lifespan -0.086519 0.1295737 -0.667717 0.5056
> Endo_ectothermendo:log.lifespan -0.234612 0.1649920 -1.421960 0.1577
>
> Correlation:
> (Intr) End_ct lg.lfs
> Endo_ectothermendo -0.580
> log.lifespan -0.798 0.676
> Endo_ectothermendo:log.lifespan 0.619 -0.912 -0.778
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -2.41100472 0.05857414 0.58071564 1.08908175 2.41709030
>
> Residual standard error: 0.7673209
> Degrees of freedom: 123 total; 119 residual## Model diagnosis
layout(matrix(c(0, 0, 0, 0,
1, 1, 2, 2,
1, 1, 2, 2,
0, 0, 0, 0), nrow = 4, ncol = 4, byrow = TRUE))
## Checking homogeneity of variance
plot(fitted(pagel.mod.red), resid(pagel.mod.red), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", type = 'n', las = 1)
grid()
par(new = TRUE)
plot(fitted(pagel.mod.red), resid(pagel.mod.red), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", las = 1)
abline(h = 0, col = "darkorange", lwd = 2)
## Checking normality
qqnorm(resid(pagel.mod.red), col = "darkgrey", type = 'n', las = 1)
grid()
par(new = TRUE)
qqnorm(resid(pagel.mod.red), col = "darkgrey", las = 1)
qqline(resid(pagel.mod.red), col = "dodgerblue", lwd = 2)
Figure 3
#png("figure3-revised.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure3-revised.pdf")
plot(log(Average_Telomere_Length_kb[Endo_ectotherm == 'ecto']) ~ log.lifespan[Endo_ectotherm == 'ecto'], data = pruned_dat, pch = 21, bg = alpha("purple", 0.5), las = 1, xlab = "Lifespan (log yrs)", ylab = "Telomere length (log kb)", type = 'n')
grid(nx = NULL, ny = NULL, col = "lightgray", lwd = 1)
par(new = TRUE)
plot(log(Average_Telomere_Length_kb[Endo_ectotherm == 'ecto']) ~ log.lifespan[Endo_ectotherm == 'ecto'], data = pruned_dat, pch = 21, bg = alpha("purple", 0.5), las = 1, xlab = "Lifespan (log yrs)", ylab = "Telomere length (log kb)")
points(log(Average_Telomere_Length_kb[Endo_ectotherm == 'endo']) ~ log.lifespan[Endo_ectotherm == 'endo'], data = pruned_dat, pch = 21, bg = alpha("orange", 0.5), las = 1, xlab = "Lifespan (log yrs)", ylab = "Telomere length (log kb)")
SSX <- sum(round((log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'ecto']) - mean(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'ecto'])))^2), 2)
s2 <- var(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'ecto']))
n <- length(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'ecto']))
x <- seq(min(pruned_dat$log.lifespan[pruned_dat$Endo_ectotherm == 'ecto']), max(pruned_dat$log.lifespan[pruned_dat$Endo_ectotherm == 'ecto']), 0.06)
m.x <- mean(round(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'ecto']), 1))
se <- sqrt(s2*((1/n) + (((x - m.x)^2)/SSX)))
is <- qt(0.975, df = n - 2)
ii <- qt(0.025, df = n - 2)
ic.s <- se*is
ic.i <- se*ii
upper.i <- (coef(mod.plot)[1] + coef(mod.plot)[2]*x) + ic.s
lower.i <- (coef(mod.plot)[1] + coef(mod.plot)[2]*x) + ic.i
##par(new = TRUE)
polygon(c(rev(x), x), c(rev(lower.i), upper.i), border = FALSE, col = alpha("purple", 0.3))
lines(x = x, y = (coef(mod.plot)[1] + (coef(mod.plot)[2]*x)), lwd = 2, col = alpha("purple", 0.5))
SSX <- sum(round((log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'endo']) - mean(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'endo'])))^2), 2)
s2 <- var(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'endo']))
n <- length(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'endo']))
x <- seq(min(pruned_dat$log.lifespan[pruned_dat$Endo_ectotherm == 'endo']), max(pruned_dat$log.lifespan[pruned_dat$Endo_ectotherm == 'endo']), 0.06)
m.x <- mean(round(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'endo']), 1))
se <- sqrt(s2*((1/n) + (((x - m.x)^2)/SSX)))
is <- qt(0.975, df = n - 2)
ii <- qt(0.025, df = n - 2)
ic.s <- se*is
ic.i <- se*ii
upper.i <- ((coef(mod.plot)[1] + coef(mod.plot)[3]) + (coef(mod.plot)[2] + coef(mod.plot)[4])*x) + ic.s
lower.i <- ((coef(mod.plot)[1] + coef(mod.plot)[3]) + (coef(mod.plot)[2] + coef(mod.plot)[4])*x) + ic.i
polygon(c(rev(x), x), c(rev(lower.i), upper.i), border = FALSE, col = alpha("orange", 0.3))
lines(x = x, y = ((coef(mod.plot)[1] + coef(mod.plot)[3]) + (coef(mod.plot)[2] + coef(mod.plot)[4])*x), lwd = 2, col = alpha("orange", 0.5))
##summary(mod.plot)
legend('bottomleft', legend = c('Ectotherms', 'Endotherms'), pch = 16, col = c('purple', 'orange'), lwd = 1, bty = 'n', cex = 0.8)
#dev.off()Figure 2
## Phylogenetic signal (Blomberg et al.'s K statitics
telo.length <- setNames(log(pruned_dat$Average_Telomere_Length_kb), rownames(pruned_dat))
telo.length> Dicentrarchus_labrax Pseudocaranx_wrighti Xiphophorus_hellerii Xiphophorus_maculatus Carassius_auratus
> 1.2354715 1.4011830 1.4469190 1.4469190 1.5260563
> Cyprinus_carpio Lutjanus_argentimaculatus Upeneichthys_vlamingii Acanthopagrus_schlegelii Macquaria_ambigua
> 1.5303947 1.6601310 1.7011051 1.7155981 1.8148247
> Nothobranchius_furzeri Nothobranchius_rachovii Gadus_morhua Platycephalus_bassensis Danio_rerio
> 1.8885837 2.0502702 2.4313300 2.5006159 2.7725887
> Alligator_sinensis Lacerta_agilis Zootoca_vivipara Emys_orbicularis Oncorhynchus_mykiss
> 2.7825391 2.9231616 2.9606231 2.9957323 2.9957323
> Neoceratodus_forsteri Urolophus_paucimaculatus Callorhinchus_milii Liasis_fuscus Meriones_unguiculatus
> 3.1871787 3.1941732 3.2023399 3.2976870 3.3672958
> Alligator_mississippiensis Chinchilla_lanigera Ondatra_zibethicus Mesocricetus_auratus Trachemys_scripta
> 3.4307562 3.6375862 3.8066625 3.9120230 3.9120230
> Myocastor_coypus Marmota_monax Tamias_striatus Accipiter_gentilis Haliaeetus_leucocephalus
> 3.9889840 4.0943446 4.1431347 0.2623643 0.7419373
> Chrysolophus_pictus Zalophus_californianus Aphelocoma_coerulescens Castor_canadensis Ateles_geoffroyi
> 1.1631508 1.6094379 1.7101878 1.9459101 1.9459101
> Fregata_minor Pygoscelis_adeliae Uria_lomvia Crocuta_crocuta Macronectes_giganteus
> 1.9586853 1.9947003 2.0412203 2.0794415 2.0893919
> Macronectes_halli Calonectris_diomedea Balaena_mysticetus Saimiri_sciureus Pteropus_rodricensis
> 2.1488510 2.1587147 2.1972246 2.1972246 2.1972246
> Homo_sapiens Aplodontia_rufa Peromyscus_maniculatus Tachycineta_albilinea Hirundo_rustica
> 2.1972246 2.1972246 2.1972246 2.2016592 2.2586332
> Haematopus_ostralegus Rissa_tridactyla Diomedea_exulans Pan_paniscus Hydrochoerus_hydrochaeris
> 2.2813615 2.2828925 2.2925348 2.3025851 2.3025851
> Giraffa_camelopardalis Pongo_pygmaeus Ceratotherium_simum Vireo_olivaceus Turdus_merula
> 2.3025851 2.3025851 2.3025851 2.3841651 2.4024304
> Larus_crassirostris Pan_troglodytes Passerculus_sandwichensis Branta_leucopsis Amazona_amazonica
> 2.4362415 2.4510051 2.4535880 2.4638532 2.4680995
> Myrmecophaga_tridactyla Tapirus_indicus Ursus_maritimus Equus_grevyi Sterna_hirundo
> 2.4849066 2.4849066 2.4849066 2.4849066 2.5494452
> Eschrichtius_robustus Chaetophractus_vellerosus Buteo_jamaicensis Aphelocoma_ultramarina Tachycineta_bicolor
> 2.5649494 2.5649494 2.5877640 2.6100698 2.6100698
> Taeniopygia_guttata Calidris_alpina Meleagris_gallopavo Loxodonta_africana Muntiacus_reevesi
> 2.6173958 2.6253930 2.6318888 2.6390573 2.6390573
> Muntiacus_muntjak Fulmarus_glacialis Cepphus_grylle Sula_sula Macaca_nemestrina
> 2.6390573 2.6567569 2.6803364 2.6946272 2.7013612
> Elephas_maximus Tupaia_tana Procavia_capensis Macaca_fascicularis Gorilla_gorilla
> 2.7080502 2.7080502 2.7080502 2.7146947 2.7343675
> Strigops_habroptila Heterocephalus_glaber Macaca_mulatta Lonchura_striata Tursiops_truncatus
> 2.7507904 2.7725887 2.7725887 2.7948393 2.8332133
> Choloepus_hoffmanni Tachymarptis_melba Lemur_catta Oceanodroma_leucorhoa Ochotona_princeps
> 2.8332133 2.9329522 2.9444390 2.9460167 2.9957323
> Riparia_riparia Colinus_virginianus Ailurus_fulgens Tadarida_brasiliensis Myotis_lucifugus
> 3.0568274 3.1527360 3.2188758 3.2580965 3.4011974
> Falco_sparverius Didelphis_virginiana Galegeeska_rufescens Macroscelides_proboscideus Atelerix_albiventris
> 3.4078419 3.5553481 3.6109179 3.6109179 3.6375862
> Buteo_swainsoni Parus_major Lepus_californicus Sciurus_carolinensis Oryctolagus_cuniculus
> 3.7471484 3.8971124 3.9120230 3.9120230 3.9120230
> Setifer_setosus Sylvilagus_aquaticus Panthera_tigris
> 3.9120230 3.9120230 3.9120230k_tl <- phylosig(pruned_data_tree, telo.length, test = TRUE, nsim = 10000)
attributes(k_tl)> $names
> [1] "K" "P" "sim.K"
>
> $class
> [1] "phylosig"
>
> $method
> [1] "K"
>
> $test
> [1] TRUE
>
> $se
> [1] FALSEhead(k_tl$sim.K)> [1] 0.16697908 0.08027691 0.10675328 0.05298757 0.05080825 0.03556077k_tl$K> [1] 0.1669791k_tl$P> [1] 2e-04#png("figure2.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure2.pdf")
plot(k_tl$sim.K, bty = "o", ylim = c(0, 3000), las = 1, ylab = "Null distribution of K", xlab = "K", main = "", type = 'n', axes = FALSE)
grid()
par(new = TRUE)
hist(k_tl$sim.K, bty = "o", ylim = c(0, 3000), las = 1, ylab = "Null distribution of K", xlab = "K", main = "")
abline(v = k_tl$K, lwd = 2, lty = "dotted")
text(x = k_tl$K - 0.02, y = 2500, "Observed value \n of K")
box()
#dev.off()
## Estimating lambda from PGLS models
pagel.model <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass, correlation = corPagel(0.8, phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(pagel.model)> Generalized least squares fit by maximum likelihood
> Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass
> Data: pruned_dat
> AIC BIC logLik
> 235.1572 252.0303 -111.5786
>
> Correlation Structure: corPagel
> Formula: ~Scientific_name
> Parameter estimate(s):
> lambda
> 0.6254488
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.629715 0.3916907 9.266788 0.0000
> Endo_ectothermendo -0.819368 0.2207192 -3.712263 0.0003
> log.lifespan -0.134538 0.1130636 -1.189932 0.2364
> log.mass -0.033764 0.0258447 -1.306413 0.1939
>
> Correlation:
> (Intr) End_ct lg.lfs
> Endo_ectothermendo -0.041
> log.lifespan -0.502 -0.101
> log.mass 0.075 0.013 -0.678
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -2.34917753 -0.02875553 0.50192730 0.92152007 2.39155284
>
> Residual standard error: 0.8117616
> Degrees of freedom: 123 total; 119 residualintervals(pagel.model, which = 'var-cov')$corStruct> lower est. upper
> lambda 0.3607961 0.6254488 0.8901015
> attr(,"label")
> [1] "Correlation structure:"pagel.model0 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name, fixed = TRUE), data = pruned_dat, method = "ML")
summary(pagel.model0)> Generalized least squares fit by maximum likelihood
> Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass
> Data: pruned_dat
> AIC BIC logLik
> 281.5872 295.6481 -135.7936
>
> Correlation Structure: corPagel
> Formula: ~Scientific_name
> Parameter estimate(s):
> lambda
> 0
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.0127181 0.27461330 10.970765 0.0000
> Endo_ectothermendo -0.0099686 0.15910482 -0.062654 0.9501
> log.lifespan -0.1146562 0.11583847 -0.989794 0.3243
> log.mass -0.0007418 0.02556906 -0.029011 0.9769
>
> Correlation:
> (Intr) End_ct lg.lfs
> Endo_ectothermendo -0.108
> log.lifespan -0.748 -0.291
> log.mass 0.131 0.098 -0.650
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -3.2551171 -0.5625105 -0.0212407 0.5130302 1.9228926
>
> Residual standard error: 0.7298431
> Degrees of freedom: 123 total; 119 residualintervals(pagel.model0, which = 'var-cov')> Approximate 95% confidence intervals
>
> Residual standard error:
> lower est. upper
> 0.6489217 0.7298431 0.8340042pagel.model0.5 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass, correlation = corPagel(0.5, phy = pruned_data_tree, form = ~Scientific_name, fixed = TRUE), data = pruned_dat, method = "ML")
summary(pagel.model0.5)> Generalized least squares fit by maximum likelihood
> Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass
> Data: pruned_dat
> AIC BIC logLik
> 233.9806 248.0415 -111.9903
>
> Correlation Structure: corPagel
> Formula: ~Scientific_name
> Parameter estimate(s):
> lambda
> 0.5
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.566355 0.3532026 10.097194 0.0000
> Endo_ectothermendo -0.752855 0.2101979 -3.581649 0.0005
> log.lifespan -0.131519 0.1105688 -1.189477 0.2366
> log.mass -0.031862 0.0254440 -1.252233 0.2129
>
> Correlation:
> (Intr) End_ct lg.lfs
> Endo_ectothermendo -0.053
> log.lifespan -0.542 -0.101
> log.mass 0.094 0.003 -0.685
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -2.57822131 -0.06855442 0.51155121 0.98428268 2.56306826
>
> Residual standard error: 0.7436376
> Degrees of freedom: 123 total; 119 residualintervals(pagel.model0.5, which = 'var-cov')> Approximate 95% confidence intervals
>
> Residual standard error:
> lower est. upper
> 0.6611868 0.7436376 0.8497675pagel.model1 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass, correlation = corPagel(1, phy = pruned_data_tree, form = ~Scientific_name, fixed = TRUE), data = pruned_dat, method = "ML")
summary(pagel.model1)> Generalized least squares fit by maximum likelihood
> Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass
> Data: pruned_dat
> AIC BIC logLik
> 272.3855 286.4464 -131.1927
>
> Correlation Structure: corPagel
> Formula: ~Scientific_name
> Parameter estimate(s):
> lambda
> 1
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) 3.901248 0.8933978 4.366753 0.0000
> Endo_ectothermendo -1.172484 0.2869553 -4.085949 0.0001
> log.lifespan -0.142480 0.1268266 -1.123421 0.2635
> log.mass -0.048054 0.0328901 -1.461046 0.1466
>
> Correlation:
> (Intr) End_ct lg.lfs
> Endo_ectothermendo -0.025
> log.lifespan -0.268 -0.136
> log.mass -0.086 0.108 -0.499
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.11576421 0.04130495 0.29665511 0.48025508 1.15040877
>
> Residual standard error: 1.926538
> Degrees of freedom: 123 total; 119 residualintervals(pagel.model1, which = 'var-cov')> Approximate 95% confidence intervals
>
> Residual standard error:
> lower est. upper
> 1.712933 1.926538 2.201488anova(pagel.model, pagel.model0)> Model df AIC BIC logLik Test L.Ratio p-value
> pagel.model 1 6 235.1572 252.0303 -111.5786
> pagel.model0 2 5 281.5872 295.6481 -135.7936 1 vs 2 48.42996 <.0001anova(pagel.model, pagel.model0.5)> Model df AIC BIC logLik Test L.Ratio p-value
> pagel.model 1 6 235.1572 252.0303 -111.5786
> pagel.model0.5 2 5 233.9806 248.0415 -111.9903 1 vs 2 0.8233961 0.3642anova(pagel.model0, pagel.model0.5)> Model df AIC BIC logLik
> pagel.model0 1 5 281.5872 295.6481 -135.7936
> pagel.model0.5 2 5 233.9806 248.0415 -111.9903anova(pagel.model, pagel.model1)> Model df AIC BIC logLik Test L.Ratio p-value
> pagel.model 1 6 235.1572 252.0303 -111.5786
> pagel.model1 2 5 272.3855 286.4464 -131.1927 1 vs 2 39.22828 <.0001pruned_dat$Telomerase_activity[pruned_dat$Telomerase_activity == 0] <- "absent"
pruned_dat$Telomerase_activity[pruned_dat$Telomerase_activity == 1] <- "present"
pruned_dat$Telomerase_activity[is.na(pruned_dat$Telomerase_activity)] <- "N/A"
tel.act <- setNames(pruned_dat$Telomerase_activity, rownames(pruned_dat))
tel.act> Dicentrarchus_labrax Pseudocaranx_wrighti Xiphophorus_hellerii Xiphophorus_maculatus Carassius_auratus
> "N/A" "N/A" "N/A" "N/A" "N/A"
> Cyprinus_carpio Lutjanus_argentimaculatus Upeneichthys_vlamingii Acanthopagrus_schlegelii Macquaria_ambigua
> "N/A" "N/A" "N/A" "N/A" "N/A"
> Nothobranchius_furzeri Nothobranchius_rachovii Gadus_morhua Platycephalus_bassensis Danio_rerio
> "present" "N/A" "present" "N/A" "present"
> Alligator_sinensis Lacerta_agilis Zootoca_vivipara Emys_orbicularis Oncorhynchus_mykiss
> "absent" "N/A" "N/A" "N/A" "present"
> Neoceratodus_forsteri Urolophus_paucimaculatus Callorhinchus_milii Liasis_fuscus Meriones_unguiculatus
> "N/A" "N/A" "N/A" "N/A" "present"
> Alligator_mississippiensis Chinchilla_lanigera Ondatra_zibethicus Mesocricetus_auratus Trachemys_scripta
> "N/A" "present" "present" "present" "N/A"
> Myocastor_coypus Marmota_monax Tamias_striatus Accipiter_gentilis Haliaeetus_leucocephalus
> "present" "present" "present" "N/A" "N/A"
> Chrysolophus_pictus Zalophus_californianus Aphelocoma_coerulescens Castor_canadensis Ateles_geoffroyi
> "N/A" "absent" "N/A" "absent" "absent"
> Fregata_minor Pygoscelis_adeliae Uria_lomvia Crocuta_crocuta Macronectes_giganteus
> "N/A" "N/A" "N/A" "absent" "N/A"
> Macronectes_halli Calonectris_diomedea Balaena_mysticetus Saimiri_sciureus Pteropus_rodricensis
> "N/A" "N/A" "absent" "absent" "absent"
> Homo_sapiens Aplodontia_rufa Peromyscus_maniculatus Tachycineta_albilinea Hirundo_rustica
> "absent" "absent" "present" "N/A" "N/A"
> Haematopus_ostralegus Rissa_tridactyla Diomedea_exulans Pan_paniscus Hydrochoerus_hydrochaeris
> "N/A" "N/A" "N/A" "absent" "absent"
> Giraffa_camelopardalis Pongo_pygmaeus Ceratotherium_simum Vireo_olivaceus Turdus_merula
> "absent" "absent" "absent" "N/A" "N/A"
> Larus_crassirostris Pan_troglodytes Passerculus_sandwichensis Branta_leucopsis Amazona_amazonica
> "N/A" "N/A" "N/A" "N/A" "N/A"
> Myrmecophaga_tridactyla Tapirus_indicus Ursus_maritimus Equus_grevyi Sterna_hirundo
> "absent" "absent" "absent" "absent" "present"
> Eschrichtius_robustus Chaetophractus_vellerosus Buteo_jamaicensis Aphelocoma_ultramarina Tachycineta_bicolor
> "absent" "absent" "N/A" "N/A" "present"
> Taeniopygia_guttata Calidris_alpina Meleagris_gallopavo Loxodonta_africana Muntiacus_reevesi
> "present" "N/A" "N/A" "absent" "absent"
> Muntiacus_muntjak Fulmarus_glacialis Cepphus_grylle Sula_sula Macaca_nemestrina
> "absent" "N/A" "N/A" "N/A" "N/A"
> Elephas_maximus Tupaia_tana Procavia_capensis Macaca_fascicularis Gorilla_gorilla
> "absent" "present" "absent" "absent" "N/A"
> Strigops_habroptila Heterocephalus_glaber Macaca_mulatta Lonchura_striata Tursiops_truncatus
> "N/A" "present" "absent" "N/A" "absent"
> Choloepus_hoffmanni Tachymarptis_melba Lemur_catta Oceanodroma_leucorhoa Ochotona_princeps
> "absent" "N/A" "absent" "present" "present"
> Riparia_riparia Colinus_virginianus Ailurus_fulgens Tadarida_brasiliensis Myotis_lucifugus
> "N/A" "N/A" "absent" "present" "present"
> Falco_sparverius Didelphis_virginiana Galegeeska_rufescens Macroscelides_proboscideus Atelerix_albiventris
> "N/A" "present" "present" "present" "present"
> Buteo_swainsoni Parus_major Lepus_californicus Sciurus_carolinensis Oryctolagus_cuniculus
> "N/A" "N/A" "absent" "present" "absent"
> Setifer_setosus Sylvilagus_aquaticus Panthera_tigris
> "present" "absent" "present"TA <- to.matrix(tel.act, unique(tel.act))
TA <- TA[pruned_data_tree$tip.label, ]
life.span <- setNames(pruned_dat$log.lifespan, rownames(pruned_dat))
log_mass <- setNames(pruned_dat$log.mass, rownames(pruned_dat))
telo.length <- setNames(log(pruned_dat$Average_Telomere_Length_kb), rownames(pruned_dat))plotTree(pruned_data_tree, ftype = "off", mar = c(3, 2, 2, 3))
tiplabels(pie = TA, piecol = c("white", "gray", "black"), cex = 0.22, offset = 4.3)
par(xpd = TRUE)
legend(x = 150, y = 0.2, legend = unique(tel.act), pch = 21, pt.bg = c("white", "gray", "black"), pt.cex = 1, bty = "n", title = "Telomerase activity", cex = 0.7, horiz = TRUE)
par(new = TRUE)
par(mar = c(3, 32, 2, 1.1))
barplot(life.span[pruned_data_tree$tip.label], horiz = TRUE, width = 1, space = 0,
ylim = c(1, length(pruned_data_tree$tip.label))-0.5, names = "", las = 2, cex.axis = 0.5, axes = FALSE)
axis(1, at = round(seq(min(life.span), max(life.span), 1.5), 1), labels = FALSE)
text(round(seq(min(life.span), max(life.span), 1.5), 1), par("usr")[3] - 0.2, labels = round(seq(min(life.span), max(life.span), 1.5), 1), srt = 50, pos = 1, xpd = TRUE, cex = 0.5, offset = 1)
mtext("Lifespan \n (years)", side = 1, line = 1.6, cex = 0.5, font = 2)
## Ancestral state reconstruction of telomere size
fit <- fastAnc(pruned_data_tree, telo.length, vars = TRUE, CI = TRUE)
fit$CI[1, ]> [1] 1.156184 4.676634obj <- contMap(pruned_data_tree, telo.length, plot = FALSE)
Figure 1
##png("figure1.png", width = 7, height = 7, units = "in", res = 360)
##pdf("figure1.pdf")
plot(obj, ftype = "off", legend = FALSE, ylim = c(1-0.09*(Ntip(obj$tree)-1), Ntip(obj$tree)), mar = c(1, 0.1, 1, 4), lwd = 1.5)
add.color.bar(150, obj$cols,title = "Log telomere length (kb)", lims = obj$lims, digits = 3, prompt = FALSE, x = 0,
y = 1-0.08*(Ntip(obj$tree)-1), lwd = 4, fsize = 0.6, subtitle = "")
par(xpd = TRUE)
text(300, 119, 'Testudines')
text(300, 101, 'Aves', col = 'red')
text(300, 82, 'Crocodilia')
text(240, 74, 'Squamata')
text(240, 40, 'Mammalia', col = 'red')
text(160, 15, 'Osteichthyes')
text(160, 4.5, 'Chondrichthyes')
#cladelabels(pruned_data_tree, node = 217, "Ectotherms", offset = 6)
#cladelabels(pruned_data_tree, node = 153, "Endotherms", offset = 6)
#segments(518, 0, 518, 18)
#segments(518, 0, 508, 0)
#segments(518, 18, 508, 18)
#text(528, 9, 'Ectotherms', srt = 90)
tiplabels(pie = TA, piecol = c("gray", "black", "white"), cex = 0.16, offset = 5.7)
legend(x = 200, y = -5, legend = unique(tel.act), pch = 21, pt.bg = c("gray", "black", "white"), pt.cex = 1, bty = "n", title = "Telomerase activity", cex = 0.7, horiz = TRUE)
par(new = TRUE)
par(mar = c(3.6, 30.8, 3, 2.2))
barplot(life.span[pruned_data_tree$tip.label], horiz = TRUE, width = 1.07, space = 0,
ylim = c(1, length(pruned_data_tree$tip.label))-0.5, names = "", las = 2, cex.axis = 0.5, axes = FALSE)
axis(1, at = round(seq(min(life.span), max(life.span), 1.5), 1), labels = FALSE)
text(round(seq(min(life.span), max(life.span), 1.5), 1), par("usr")[3] - 0.2, labels = round(seq(min(life.span), max(life.span), 1.5), 1), srt = 50, pos = 1, xpd = TRUE, cex = 0.5, offset = 1)
mtext("Log lifespan \n (years)", side = 1, line = 1.6, cex = 0.5, font = 2)
##dev.off()## Correlated evolution under the threshold model
## Removing NAs from the dataset
pruned_dat$Telomerase_activity> [1] "N/A" "N/A" "N/A" "N/A" "N/A" "N/A" "N/A" "N/A" "N/A" "N/A" "present" "N/A" "present"
> [14] "N/A" "present" "absent" "N/A" "N/A" "N/A" "present" "N/A" "N/A" "N/A" "N/A" "present" "N/A"
> [27] "present" "present" "present" "N/A" "present" "present" "present" "N/A" "N/A" "N/A" "absent" "N/A" "absent"
> [40] "absent" "N/A" "N/A" "N/A" "absent" "N/A" "N/A" "N/A" "absent" "absent" "absent" "absent" "absent"
> [53] "present" "N/A" "N/A" "N/A" "N/A" "N/A" "absent" "absent" "absent" "absent" "absent" "N/A" "N/A"
> [66] "N/A" "N/A" "N/A" "N/A" "N/A" "absent" "absent" "absent" "absent" "present" "absent" "absent" "N/A"
> [79] "N/A" "present" "present" "N/A" "N/A" "absent" "absent" "absent" "N/A" "N/A" "N/A" "N/A" "absent"
> [92] "present" "absent" "absent" "N/A" "N/A" "present" "absent" "N/A" "absent" "absent" "N/A" "absent" "present"
> [105] "present" "N/A" "N/A" "absent" "present" "present" "N/A" "present" "present" "present" "present" "N/A" "N/A"
> [118] "absent" "present" "absent" "present" "absent" "present"pruned_dat_not_NAs <- pruned_dat[!pruned_dat$Telomerase_activity == "N/A", ]
pruned_dat_not_NAs$Telomerase_activity> [1] "present" "present" "present" "absent" "present" "present" "present" "present" "present" "present" "present" "present" "absent"
> [14] "absent" "absent" "absent" "absent" "absent" "absent" "absent" "absent" "present" "absent" "absent" "absent" "absent"
> [27] "absent" "absent" "absent" "absent" "absent" "present" "absent" "absent" "present" "present" "absent" "absent" "absent"
> [40] "absent" "present" "absent" "absent" "present" "absent" "absent" "absent" "absent" "present" "present" "absent" "present"
> [53] "present" "present" "present" "present" "present" "absent" "present" "absent" "present" "absent" "present"check2 <- name.check(pruned_data_tree, pruned_dat_not_NAs)
rm_phy2 <- check2$tree_not_data
rm_dat2 <- check2$data_not_tree
pruned_data_tree_not_NAs <- drop.tip(pruned_data_tree, rm_phy2)
pruned_dat_not_NAs <- subset(pruned_dat_not_NAs, subset = pruned_dat_not_NAs$Scientific_name %in% pruned_data_tree_not_NAs$tip, select = names(pruned_dat_not_NAs))
pruned_dat_not_NAs$Telomerase_activity <- as.factor(pruned_dat_not_NAs$Telomerase_activity)
str(pruned_dat_not_NAs)> 'data.frame': 63 obs. of 29 variables:
> $ Common.name : chr "African killifish" "Atlantic cod" "Zebrafish" "Chinese alligator" ...
> $ Scientific_name : chr "Nothobranchius_furzeri" "Gadus_morhua" "Danio_rerio" "Alligator_sinensis" ...
> $ Domain : chr "Eukaryota" "Eukaryota" "Eukaryota" "Eukaryota" ...
> $ Kingdom : chr "Metazoa" "Metazoa" "Metazoa" "Metazoa" ...
> $ Phylum : chr "Chordata" "Chordata" "Chordata" "Chordata" ...
> $ Class : chr "Fish" "Fish" "Fish" "Reptilia" ...
> $ Order : chr "Cyprinodontiformes" "Gadiformes" "Cypriniformes" "Crocodilia" ...
> $ Family : chr "Nothobranchiidae" "Gadidae" "Cyprinidae" "Alligatoridae" ...
> $ Genus : chr "Nothobranchius" "Gadus" "Danio" "Alligator" ...
> $ Species : chr "furzeri" "morhua" "rerio" "sinensis" ...
> $ Endo_ectotherm : chr "ecto" "ecto" "ecto" "ecto" ...
> $ Adult_mass_grams : num 2 52800 0.5 14600 12235 ...
> $ log_mass : num 0.301 4.723 -0.301 4.164 4.088 ...
> $ Lifespan_years : num 0.23 25 5.5 65 11 6.3 17.2 10 3.9 19 ...
> $ Loglifespan : num -0.638 1.398 0.74 1.813 1.041 ...
> $ Average_Telomere_Length_kb: num 6.61 11.37 16 16.16 20 ...
> $ Telomerase_activity : Factor w/ 2 levels "absent","present": 2 2 2 1 2 2 2 2 2 2 ...
> $ Tissue.type.for.TA : chr "Various tissues" "Multiple somatic tissues" "Various tissues" "Blood" ...
> $ Source.for.TA : chr "Hartmann et al. 2009" "Abechuco et al. 2014" "Achelin et al. 2011" "Guo et al. 2023" ...
> $ Tissue.type.for.TL : chr "Muscle and skin" "Multiple somatic tissues" "Various tissues" "Blood" ...
> $ Methodology.for.TL : chr "TRF (denatured)" "TRF (denatured)" "TRF (denatured)" "TRF (denatured)" ...
> $ TRF_type : chr "A" "A" "A" "A" ...
> $ Source.for.TL : chr "Hartmann et al. 2009" "de Abechuco et al. 2016" "Lund et al. 2009; Achelin et al. 2011" "Xu et al. 2009" ...
> $ Source.for.Lifespan : chr "Valdesalici and Cellerino 2003" "AnAge" "AnAge" "AnAge" ...
> $ Source.for.Mass : chr "Pola?ik and Reichard 2010" "AnAge" "Fishbase" "AnAge" ...
> $ Captivity : int 0 0 0 0 0 0 0 0 0 0 ...
> $ Notes : chr "" "" "" "" ...
> $ log.lifespan : num 0.207 3.258 1.872 4.19 2.485 ...
> $ log.mass : num 1.099 10.874 0.405 9.589 9.412 ...head(pruned_dat_not_NAs)> Common.name Scientific_name Domain Kingdom Phylum Class Order Family
> Nothobranchius_furzeri African killifish Nothobranchius_furzeri Eukaryota Metazoa Chordata Fish Cyprinodontiformes Nothobranchiidae
> Gadus_morhua Atlantic cod Gadus_morhua Eukaryota Metazoa Chordata Fish Gadiformes Gadidae
> Danio_rerio Zebrafish Danio_rerio Eukaryota Metazoa Chordata Fish Cypriniformes Cyprinidae
> Alligator_sinensis Chinese alligator Alligator_sinensis Eukaryota Metazoa Chordata Reptilia Crocodilia Alligatoridae
> Oncorhynchus_mykiss Rainbow trout Oncorhynchus_mykiss Eukaryota Metazoa Chordata Fish Salmoniformes Salmonidae
> Meriones_unguiculatus Mongolian Gerbil Meriones_unguiculatus Eukaryota Metazoa Chordata Mammalia Rodentia Muridae
> Genus Species Endo_ectotherm Adult_mass_grams log_mass Lifespan_years Loglifespan
> Nothobranchius_furzeri Nothobranchius furzeri ecto 2.0 0.301030 0.23 -0.6382722
> Gadus_morhua Gadus morhua ecto 52800.0 4.722634 25.00 1.3979400
> Danio_rerio Danio rerio ecto 0.5 -0.301030 5.50 0.7403627
> Alligator_sinensis Alligator sinensis ecto 14600.0 4.164353 65.00 1.8129134
> Oncorhynchus_mykiss Oncorhynchus mykiss ecto 12235.0 4.087604 11.00 1.0413927
> Meriones_unguiculatus Meriones unguiculatus ecto 53.2 1.725912 6.30 0.7993405
> Average_Telomere_Length_kb Telomerase_activity Tissue.type.for.TA Source.for.TA
> Nothobranchius_furzeri 6.610 present Various tissues Hartmann et al. 2009
> Gadus_morhua 11.374 present Multiple somatic tissues Abechuco et al. 2014
> Danio_rerio 16.000 present Various tissues Achelin et al. 2011
> Alligator_sinensis 16.160 absent Blood Guo et al. 2023
> Oncorhynchus_mykiss 20.000 present Multiple tissues Panasiak et al. 2023
> Meriones_unguiculatus 29.000 present Multiple somatic tissues Seluanov et al. 2007
> Tissue.type.for.TL Methodology.for.TL TRF_type Source.for.TL
> Nothobranchius_furzeri Muscle and skin TRF (denatured) A Hartmann et al. 2009
> Gadus_morhua Multiple somatic tissues TRF (denatured) A de Abechuco et al. 2016
> Danio_rerio Various tissues TRF (denatured) A Lund et al. 2009; Achelin et al. 2011
> Alligator_sinensis Blood TRF (denatured) A Xu et al. 2009
> Oncorhynchus_mykiss Blood TRF (denatured) A Lejnine et al. 1995
> Meriones_unguiculatus Multiple somatic tissues TRF (denatured) A Seluanov et al. 2007
> Source.for.Lifespan Source.for.Mass Captivity Notes log.lifespan log.mass
> Nothobranchius_furzeri Valdesalici and Cellerino 2003 Pola?ik and Reichard 2010 0 0.2070142 1.0986123
> Gadus_morhua AnAge AnAge 0 3.2580965 10.8742854
> Danio_rerio AnAge Fishbase 0 1.8718022 0.4054651
> Alligator_sinensis AnAge AnAge 0 4.1896547 9.5888453
> Oncorhynchus_mykiss AnAge AnAge 0 2.4849066 9.4121377
> Meriones_unguiculatus AnAge AnAge 0 1.9878743 3.9926809pruned_data_tree_not_NAs>
> Phylogenetic tree with 63 tips and 62 internal nodes.
>
> Tip labels:
> Oncorhynchus_mykiss, Nothobranchius_furzeri, Gadus_morhua, Danio_rerio, Didelphis_virginiana, Atelerix_albiventris, ...
> Node labels:
> '12', '48', '57', '49', '194', '166', ...
>
> Rooted; includes branch lengths.name.check(pruned_data_tree_not_NAs, pruned_dat_not_NAs)> [1] "OK"names(pruned_dat_not_NAs)> [1] "Common.name" "Scientific_name" "Domain" "Kingdom"
> [5] "Phylum" "Class" "Order" "Family"
> [9] "Genus" "Species" "Endo_ectotherm" "Adult_mass_grams"
> [13] "log_mass" "Lifespan_years" "Loglifespan" "Average_Telomere_Length_kb"
> [17] "Telomerase_activity" "Tissue.type.for.TA" "Source.for.TA" "Tissue.type.for.TL"
> [21] "Methodology.for.TL" "TRF_type" "Source.for.TL" "Source.for.Lifespan"
> [25] "Source.for.Mass" "Captivity" "Notes" "log.lifespan"
> [29] "log.mass"## Set the number of generations
ngen <- 5e6
## Run MCMC
pruned_dat_not_NAs$Average_Telomere_Length_kb <- log(pruned_dat_not_NAs$Average_Telomere_Length_kb)
mcmc.model <- threshBayes(pruned_data_tree_not_NAs, pruned_dat_not_NAs[ , c(16, 17)], type = c("cont", "disc"), ngen = ngen,
plot = FALSE, control = list(print.interval = 5e+05))> Starting MCMC....
> generation: 500000; mean acceptance rate: 0.23
> generation: 1000000; mean acceptance rate: 0.24
> generation: 1500000; mean acceptance rate: 0.27
> generation: 2000000; mean acceptance rate: 0.2
> generation: 2500000; mean acceptance rate: 0.26
> generation: 3000000; mean acceptance rate: 0.23
> generation: 3500000; mean acceptance rate: 0.21
> generation: 4000000; mean acceptance rate: 0.23
> generation: 4500000; mean acceptance rate: 0.21
> generation: 5000000; mean acceptance rate: 0.22
> Done MCMC.mcmc.model>
> Object of class "threshBayes" consisting of a matrix (L) of
> sampled liabilities for the tips of the tree & a second matrix
> (par) with the sample model parameters & correlation.
>
> Mean correlation (r) from the posterior sample is: 0.61357.
>
> Ordination of discrete traits:
>
> Trait 2: absent <-> present## Pull out the post burn-in sample and compute HPD
r.mcmc <- tail(mcmc.model$par$r, 0.8*nrow(mcmc.model$par))
class(r.mcmc) <- "mcmc"
hpd.r <- HPDinterval(r.mcmc)
hpd.r> lower upper
> var1 0.3487216 0.8541938
> attr(,"Probability")
> [1] 0.9500012
Figure S3
## Profile plots from a Bayesian MCMC analysis of the threshold model
#png("figureS3.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figureS3.pdf")
layout(matrix(c(0, 0, 0, 0,
1, 1, 1, 1,
1, 1, 1, 1,
2, 2, 2, 2,
2, 2, 2, 2,
0, 0, 0, 0), nrow = 6, ncol = 4, byrow = TRUE))
plot(mcmc.model$par$logL ~ mcmc.model$par$gen, col = 'gray', type = 'n', ylab = 'Log(L)', xlab = 'Generation', las = 1)
grid()
par(new = TRUE)
plot(mcmc.model$par$logL ~ mcmc.model$par$gen, col = 'gray', type = 'l', ylab = 'Log(L)', xlab = 'Generation', las = 1)
mtext('(a)', side = 2, line = 1.5, at = -160, las = 1)
plot(mcmc.model$par$r ~ mcmc.model$par$gen, col = 'purple', type = 'n', ylab = 'r', xlab = 'Generation', las = 1)
grid()
par(new = TRUE)
plot(mcmc.model$par$r ~ mcmc.model$par$gen, col = 'purple', type = 'l', ylab = 'r', xlab = 'Generation', las = 1)
mtext('(b)', side = 2, line = 1.5, at = 1.2, las = 1)
#dev.off()
Figure 4
## Plot posterior density
#png("figure4.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure4.pdf")
par(las = 1)
plot(NA, xlim = c(-1, 1.5), ylim = c(0, 3), axes = FALSE, xlab = '', ylab = '')
grid()
par(new = TRUE)
plot(density(mcmc.model), xlim = c(-1, 1.5))
box()
## add whiskers to show HPD
h <- 0-par()$usr[3]
lines(x = hpd.r, y = rep(-h/2, 2))
lines(x = rep(hpd.r[1], 2), y = c(-0.3, -0.7)*h)
lines(x = rep(hpd.r[2], 2), y = c(-0.3, -0.7)*h)
#dev.off()
Sensitivity analysis - Added by Derek
Benson
## Create sensitivity analysis object using 'samp' method
samp <- samp_phylm(Average_Telomere_Length_kb ~ log.lifespan, phy = full_data_tree, data = pruned_dat, n.sim = 1000, track = FALSE)
summary(samp)> % Species Removed % Significant Intercepts Mean Intercept Change (%) Mean sDIFintercept % Significant Estimates
> 1 10 99.9 3.3375 -0.09634809 99.8
> 2 20 99.4 5.4386 -0.20865625 98.9
> 3 30 99.1 7.2334 -0.32696877 97.3
> 4 40 97.9 9.3075 -0.42784789 91.1
> 5 50 98.3 11.4862 -0.49835594 81.3
> Mean Estimate Change (%) Mean sDIFestimate
> 1 6.0234 0.011289138
> 2 9.9227 0.027871516
> 3 12.9821 0.052113105
> 4 16.9222 0.099651154
> 5 20.9871 0.006694003sensi_plot(samp, graphs = 1, param = "estimate")
sensi_plot(samp, graphs = 2, param = "estimate")
sensi_plot(samp, graphs = 3, param = "estimate")
sensi_plot(samp, graphs = 4, param = "estimate")